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Abstract 

A low-dimensional model of large Prandtl-number (P) Rayleigh Benard convection is constructed 
using some of the important modes of pseudospectral direct numerical simulations. A detailed 
bifurcation analysis of the low-dimensional model for P = 6.8 and aspect ratio of 2\/2 reveals a 
rich instability and chaos picture: steady rolls, time-periodicity, quasiperiodicity, phase locking, 
chaos, and crisis. Bifurcation analysis also reveals multiple co-existing attractors, and a window 
with time-periodicity after chaos. The results of the low-dimensional model matches quite closely 
with some of the past simulations and experimental results where they observe chaos in RBC 
through quasiperiodicity and phase locking. 
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I. INTRODUCTION 



Rayleigh Benard convection (RBC) exhibits a host of complex phenomena. Instabilities, 
patterns, and chaos are observed near the onset of convection, while turbulence is seen far 
beyond the onset Bifurcation diagrams are often used to study instabilities and chaos. 
In this paper we report a bifurcation analysis for large-Prandtl number RBC flow obtained 
using a low-dimensional model. 

The two most important parameters for RBC are the Rayleigh number R (ratio of buoy- 
ancy and dissipative force) and the Prandtl number P (ratio of viscosity and thermal dif- 
fusivity). Convection starts at the critical Rayleigh number R c , which is independent of P. 
For non-zero Prandtl number, the primary instability at the onset of convection is always 
in the form of two-dimensional (2D) straight rolls [2]. These rolls become unstable through 
secondary instabilities, and they bifurcates into a sequence of dynamic patterns {3, 4, 3|. 
Some examples of the resulting patterns are squares and hexagons, asymmetric patterns, 
oscillating patterns, relaxation oscillations of rolls and squares, etc. 6|. 

The sequence of instabilities and onset of chaos are quite different for low-Prandtl num- 
ber (low-P) convection [7J and large-Prandtl number (large-P) convection j^, t], 10, [ll]. For 



low-P convection, the 2D rolls become unstable close to the onset and wavy rolls are gen- 
erated through secondary instabilities, thus making the flow three-dimensional (3D). These 
bifurcations and route to chaos for low-P and zero-Prandtl (zero-P) number convection have 
been studied extensively (see Q, 13] and references therein). The scenario however is 
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quite different for large-P convection. The secondary instabilities are delayed here, and the 
2D rolls continue to be solutions till larger Rayleigh numbers. It has been reported tha t 2D 
convection results have significant similarities with 3D results for large-P convection 
We exploit this observation to analyze bifurcation and chaos for large-P convection using a 
low-dimensional model containing only 2D modes. A major advantage of this simplification 
is that the number of modes required for 2D convection is much fewer than 3D convection, 
thus enabling the bifurcation analysis. 

n 

Krishnamurti [16j performed extensive convection experiments on mercury (P « 0.02), 
air (P « 0.7), water (P « 6.8), freon 113 (P « 7), and silicon oil (P ~ 100). She studied 
transition from 2D convection to 3D convection and subsequent generation of oscillatory, 
chaotic, and turbulent convection. Busse and Whitehead s| also reported 'zigzag instability' 
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and 'cross-roll instability' in an experiment on silicon oil. 

Libchaber et al. [7] studied the routes to chaos for RBC in mercury (a low-P fluid) as a 
function of R and the applied mean magnetic field; at different values of the parameters they 
observed chaos through various routes: period doubling, quasiperiodic, and soft instability. 
Gollub and Benson [lOj studied route to chaos in RBC of water at two different Prandtl 
numbers, P = 2.5 and 5, for different aspect ratios T and observed very rich behaviour. For 
P = 5 and T = 3.5 they observed steady rolls till r = R/ R c = 27.2 (r is called 'reduced 
Rayleigh number'), at which point periodic flow starts. At r = 32, a second frequency 
appears in the system and the flow becomes quasiperiodic. Phase locking occurs at r = 44.4 
that finally leads to chaos at r = 46.0. Gollub and Benson also observed period doubling 
route to chaos for P = 2.5 and T = 3.5, and quasiperiodic route to chaos with the third 
frequency for P = 5 and T = 2.4. They also observed intermittency in their system. 

In another experiment, Maurer and Libchaber 17J observed frequency locking and subse- 
quent generation of chaos in liquid helium as a result of generation of new frequency modes. 
Giglio et al. HI observed period-doubling route to chaos in their convective experiment on 

r n n 

water. Berge et al. [9| found intermittency in RBC of silicon oil. Ciliberto and Rubio [18( 
reported localized oscillations and travelling waves in RBC. Morris et al. [19| discovered 
spatio-temporal chaos in RBC of silicon oil. These results indicate complex nonlinear dy- 
namics including chaos in RBC. 

Direct numerical simulation (DNS) of 2D and 3D RBC have been used to study various 
convective states including turbulence. Curry et al. [20j performed detailed 3D DNS for 
P = 10 under free-slip boundary conditions; they reported steady convection till r pa 40, 
after which single frequency oscillations are observed till r pa 45. Subsequently they ob- 
served quasiperiodicity (r pa 45 — 55), phase locking (r pa 55 — 65), and chaos (r > 65). 



Yahata 



2 if ] performed DNS using finite difference scheme (MAC method) on no-slip bound- 



ary condition for P = 5, T x = 3.5 and T y = 2. The results show a series of bifurcations 
from monoperiodic — >• biperiodic — > frequency-locked state — > chaotic state. Mukutmoni 
and Yang 22| reported a numerical study of RBC in a rectangular enclosure with insulated 
sidewalls. They observed a period-2 response after a periodic solution but the route to chaos 
is through quasiperiodicity. However, on imposing symmetry of the velocity and tempera- 
ture field about the mid-planes, they observed a period-doubling route to chaos. They have 
also reported periodic solutions after chaos. 
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In two dimensions, Moore and Weiss 



231 ] simulated P = 6.8 RBC using a spectral method 



with free-slip boundary conditions; they studied heat transport as a function of the Rayleigh 
number. McLaughlin and Orszag 24( considered RBC in air (P = 0.71) with no-slip bound- 
ary conditions; they obtained periodic, quasiperiodic, and chaotic states for Rayleigh num- 
bers between 6500 and 25000. Curry et al. 20j performed detailed DNS for P = 6.8 fluid 
in 2D and observed oscillations with a single frequency at r pa 50, and with two frequencies 
at r pa 290. They also reported weak chaos beyond r = 290, and a periodic solution after 
r pa 800. Goldhirsch et al. 



25| also simulated 2D RBC and observed complex behaviour. 



Recently, Paul et al. 26] performed 2D DNS for free-slip boundary condition for a large 
range of Rayleigh numbers and obtained steady convection (r = — 80), periodicity (r = 
80 — 660), quasiperiodicity (r = 660 — 770), and chaos (r = 770 — 890). They also observed 



periodic and steady convection beyond the chaotic state. Paul et al. [26j's results are in 
general agreement with those of Curry et al. 20] and Goldhirsch et al. with some 
difference in the Rayleigh numbers. One noticeable difference between 2D RBC and 3D 
RBC is that the secondary instabilities in 2D RBC occur at significantly higher values of r 
than those in 3D RBC. 

The origin of various convective patterns and chaos in DNS is not apparent due to various 



interactions among the large number o 
dimensional models of RBC. Curry 



modes. Rather, they have been analyzed using low- 
271 ] constructed a 14-mode model of RBC with a small 
amplitude periodic modulations in the heat equation. He observed chaos for P = 10. The 
low-dimensional model of Curry shows features similar to the experiments of Gollub and 
Benson [10] namely periodicity, quasiperiodicity, and chaos. Curry 28] also studied the 14- 
mode model without any modulation, and compared its results with those from the Lorenz 
model. 

Yahata 29] studied transition to chaos in RBC using a 48-mode system of equations 



under no-slip boundary conditions. For P = 5, T x = 2 and Y y = 3.5, he obtained periodic 
— > quasiperiodic motion with two fundamental frequencies — > quasiperiodicity with three 
frequencies — > chaos. The whole sequence of bifurcations occur in the range of r = 39.77 to 



41.04. Yahata 



30] continued the above analysis for P = 2.5 with the same aspect ratio and 



reported period-doubling route to chaos for r in the range of 24.46 to 29.35. 



In the above work, Curry 



271 ] and Yahata 29_|, |30j numerically integrate the low- 



dimensional models for certain r values and observe various patterns. Some of the short- 
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comings of this method are that it could miss the behaviour in narrow windows, and the 
bifurcation points cannot be precisely located. In the present paper we numerically time 
advance the fixed points, limit cycles, and chaotic attractors, which provides a much more 
detailed bifurcation picture than that obtained by studying patterns at selected r values. 
The approach in the present paper is very similar to a recent work by Pal et al. 13] where a 
detailed bifurcation diagram was constructed for zero-P convection using a 13-mode system; 
this system was obtained using the energetic modes of DNS. Note that in large-P convection 
chaos sets in at a much larger Rayleigh number than in low-P and zero-P convection. As a 
consequence, large-P convective flows contain a large number of energetic modes at the onset 
of chaos. A bifurcation analysis of these large number of modes is very difficult and imprac- 
tical. To circumvent this difficulty we study 2D convection, which is not a serious limitation 
for large-P convection since 2D and 3D convection have significant similarities HQ. 

We perform our low- dimensional analysis for P = 6.8, which is the Prandtl number of 
water, a representative of large-P fluid. A significant advantage of choosing this Prandtl 
number is that a large number of DNS and experiments have been performed for water. We 
find that optimal number of real Fourier modes of our low-dimensional model is 30. The 
convective patterns-steady, periodic, quasiperiodic, and chaotic rolls-observed in DNS and 
experiments are captured in our low-dimensional model. 

The outline of the paper is as follows. In Section II we describe the governing equations 
and the low- dimensional model. Details of the various convective regimes, bifurcation sce- 
nario and route to chaos associated with this model is presented in section III. The last 
section contains discussions and conclusions. 



II. THE LOW-DIMENSIONAL MODEL AND ITS NUMERICAL SIMULATION 

In RBC, a thin layer of fluid is confined between two thermally conducting horizontal 
plates that are separated by a distance d. The fluid has kinematic viscosity u, thermal 
diffusivity k, and thermal expansion coefficient a. An adverse temperature gradient (3 = 
AT/d is imposed across the fluid layer, where AT is the temperature difference across the 
layer. We assume Boussinesq approximation for the fluid [l]. The relevant hydrodynamic 
equations are nondimensionalized using the length scale d, the large-scale velocity scale 
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a/ a(3gd 2 , and the temperature scale AT to yield [l| 



/ P 

d t w + (v • V)v = -Vp + 9z + \ ^V 2 v, (1) 

V ti 

dt 9 + (v • V)6 = v 3 + ^V 2 #, (2) 
V-v = 0, (3) 

where v = (1)1,1)2,1)3) is the velocity fluctuation, 9 is the perturbations in the temperature 
field from the steady conduction state, R = ag/3d 4 /i>K is the Rayleigh number, P = u/k is 
the Prandtl number, g is the acceleration due to gravity, and z is the buoyancy direction. 

The above equations are often solved using direct numerical simulation (DNS). One of 
the popular numerical technique is pseudospectral method in which the velocity and the 
temperature fields are expanded in the Fourier/Chebyshev basis. These numerical simu- 
lations have been able to reproduce various patterns, chaos, and turbulence observed in 
experiments. The convection simulations are however very expensive in terms of computer 
time and memory. Also a large number of modes present in DNS obscures the internal dy- 
namics. A popular method to analyze such systems is a bifurcation analysis of appropriate 
low-dimensional systems. Using this technique we can study the origin of various patterns 
and chaos in RBC For our low- dimensional model we choose fourteen complex modes and 
two real modes that represent the large-scale flow structures. Expansion of the vertical 
velocity field v 3 and the temperature field 9 using these modes yields 

v 3 (x, z, t) = W W i(t) exp(ik c x) sin(7T2;) + W W3 (t) exp(ik c x) sin(37rz) 

+ Wios(t) exp(ik c x) sin(57rz) + Wwitit) exp(2ik c x) sin(27rz) 

+ W 3 Qi(t) exp(3ik c x) sm(nz) + W 3 o 3 (t) exp(3ik c x) sin(37r^) 

+ W5oi(t) exp(5ik c x) sm(irz) + c.c, 

v 2 (x,z,t) = 0, 

9(x,z,t) = 9 W i(t) exp(ik c x) sm(rrz) + 9 W3 (t) exp(ik c x) sm(37iz) 

+ $105 (t) exp(ik c x) sin(57T2;) + #202 (t) exp(2ik c x) sin(27T2;) 

+ $301 (t) exp(3ik c x) sin(7r,2) + 9 303 (t) exp(3ik c x) sin(37r£) 

+ $501 (t) exp(5ik c x) sin(7rz) + c.c. 

+ 0QO2 (t) sin(27T2;) + 6<m(t) sin(47r,2) (4) 



6 



where c.c. stands for the complex conjugate, and the three subscripts denote the Fourier 
wavenumber indices along x, y, and z directions respectively. These modes correspond to 
the free-slip boundary condition for the velocity modes. Note that v%(x, z) can be computed 
using the incompressibility condition V ■ v = 0. We choose k c = hence the aspect 

ratio of our model is 2y2. The Galerkin projection of Eqs. ([1113]) on these modes yields a set 
of thirty coupled ordinary differential equations (ODEs) for the real and imaginary parts of 
the Fourier modes. These thirty nonlinear ODEs comprise our low-dimensional model. 
The above modes represent two-dimensional rolls. It has been reported earlier that the 



2D and 3D convection have significant similarity for large-P flows [15J]. Therefore we expect 
our low-dimensional model to capture the dynamics of large-P convection. Our model with 
2D rolls has 30 modes while a full 3D low-dimensional model would have many more modes 
that would make the bifurcation analysis of the model very difficult. Note that three- 
dimensional patterns like squares are not accessible to our model. However quasiperiodicity 
and the origin of chaos are expected to be common for both 2D and 3D convection for large 
Prandtl number flows. 

[ — I 

Our low- dimensional model is quite similar to that of Curry [271 ] . A major difference is 
that in our model all the modes except #002 and #004 are complex in contrast to Curry's 
model in which they are all real. Also, we keep the mode (105), whereas Curry keeps (204). 



Several of the patterns and chaos reported in the experiments of Gollub and Benson [10] 
have been observed by Curry when he includes small amplitude modulation. We do not 
require any modulation or any additional forcing (other than buoyancy) in our model to 



produce these patterns and chaos. Yahata [23, [30|] studied RBC under no-slip boundary 



condition by expanding the velocity and temperature fields using mixed basis functions 
(Chebyshev along the buoyancy direction and Fourier along the horizontal directions) and 
observed similar behaviour. Surprisingly the patterns and chaos reported for the no-slip and 
the free-slip boundary conditions are quite similar. 

We numerically solve the low-dimensional model using random initial conditions. In 
our low-dimensional model, we observe various patterns: steady convection, periodicity, 
quasiperiodicity, and chaos at different values of Rayleigh numbers (see Fig. [1]). Figured]) 
also shows that the system becomes periodic after chaos, and then it becomes chaotic again. 

A curious feature of our model is that for the periodic solutions (r = 27.6 — 40.3), the 
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FIG. 1: Time series of the amplitude of the complex mode Wioi generated by the low-dimensional 
model at various representative reduced Rayleigh numbers (r). We observe (a) steady convection 
(r = 27.5), (b) time-periodic convection (r = 35), (c) quasi-periodicity (r = 42) and (d) chaos 
(r = 45). Subsequently, (e) a window of time-periodic state (r = 47.2) followed by (f) a chaotic 
state (r = 49) is observed. 

frequency of the Fourier amplitude is twice that of its phase (see Fig. [2]). This feature can 
also be understood using the bifurcation analysis that will be discussed below. 

The origin of the observed patterns can be studied more rigorously using the bifurcation 
analysis which is the subject of the next section. 

III. BIFURCATION ANALYSIS OF THE LOW-DIMENSIONAL MODEL 

In Fig. [3] we present a bifurcation diagram obtained by numerical integration of the low- 
dimensional model for P = 6.8 and the aspect ratio of 2y/2. The reduced Rayleigh number 
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FIG. 2: Time-series of amplitude and phase of the complex Wioi mode at r = 30. The amplitude 
oscillates at double the frequency of its phase. 

r is the bifurcation parameter in our analysis. The unstable solutions have not been shown 
in Fig. [3l To generate this bifurcation diagram, numerical simulations have been performed 
with a fixed initial condition till t = 20000 (large-scale eddy turnover time). Transients 
till t = 5000 are eliminated and the extremum values of |Wioi| are plotted for later time. 
The stability and bifurcations of the steady states in this numerically generated bifurcation 
diagram are complimented by an eigenvalue analysis of the jacobian evaluated at the fixed 
points, and the eigenvalue of the associated Floquet matrix for the limit cycles. For this 
complimentary analysis, a fixed point is obtained numerically using the Newton-Raphson 
method for a given r, and the branches of the fixed points are subsequently obtained using 
a fixed arc- length based continuation scheme (similar to the analysis in Pal et al. [3]). 

For r < 1, there is no convection in the system and heat is transported solely by conduc- 
tion. The conduction state corresponds to the trivial fixed point of our system. At r = 1, 
the system undergoes a pitch-fork bifurcation and time-independent convection states are 
born as non-trivial fixed points. Note that there are two nonzero solutions near the onset of 
convection (e.g, iM^ioi). In Fig. [3] we plot |Wioi| vs. r. The new stable roll solutions (blue 
curve labeled 'FP') remain stable till r « 27.6. 



9 




FIG. 3: Bifurcation diagram of the low-dimensional model representing large-Prandtl number RBC. 
'FP' (blue curve) is the steady roll, 'OS' (red curve) is the time-periodic roll, 'QP' (green patch) 
is the quasi-periodic roll, and 'CH' (pink patch) is the chaotic state. 'NS' indicates the Neimark- 
Sacker bifurcation point. A window of periodic and quasiperiodic states is observed in the band of 
r = 46.2 - 48.4. 

The branch of fixed points corresponding to the steady convective rolls undergoes a 
supercritical Hopf bifurcation at r = 27.6. As a consequence, the time-independent steady 
state solution becomes unstable and a new stable time-periodic state (limit cycle) is born. 
The limit cycle solution is shown as a red line with the label 'OS' in the bifurcation diagram 
(Fig. [3]). The two lines of 'OS' state designate the maxima and minima of |Wioi| respectively. 
The time variation of the modes is however more complex. As shown in Fig. [2J the amplitude 
of the mode Wioi vary with frequency twice that of its phase. This phenomena can be 
understood as follows. At the Hopf bifurcation point, the eigenvectors associated with 
the pair of purely imaginary eigenvalues ±iu have components only along the imaginary 
part of the Fourier modes (e.g., Q^Wioi))- Hence, Qf(Wioi) oscillates with the frequency uj 
corresponding to the Hopf point. The real parts of the Fourier modes are generated purely 
due to the quadratic nonlinearities involving products of two imaginary parts of the modes; 
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hence 2a; (superharmonic) is the leading frequency of the real parts and the amplitudes of 
modes. 

We determine the stability of the above time-periodic state using the Floquet theory. We 
numerically construct the fundamental (Floquet) matrix associated with the time-periodic 
state and compute its eigenvalues (called 'Floquet multipliers'). All the Floquet multipliers 
for a stable limit cycle have magnitude less than one. For an unstable limit cycle, at least 
one of them has a magnitude greater than one. When the Floquet multipliers cross along the 
positive real axis, new limit cycles may appear or disappear ('pitchfork' or 'turning point' 
bifurcation). If they cross the negative real axis, the frequency of the limit cycle doubles in 
a 'period-doubling' bifurcation. However, a pair of complex-conjugate multipliers may also 
cross the unit circle in the complex plane, wherein another frequency is generated. This 
bifurcation is known as 'Niemark-S acker' (NS) or a secondary Hopf bifurcation. 

Fig. H] illustrates the magnitude of the largest Floquet multiplier as a function of the 
reduced Rayleigh number r, while Fig. [5] shows the movement of the Floquet multipliers 
around several values of r. For these calculations, we proceed as follows. The limit cycle 
along with its time-period is obtained as a fixed point of an appropriately defined map 
(described in appendix [A} using the Newton- Raphson method for a given r. The branch of 
limit cycles is subsequently computed using a fixed arc-length based continuation scheme. 
The fundamental matrix and the eigenvalues associated with the limit cycles for each value 
of r are then evaluated numerically. 

Our computations reveal that the largest Floquet multiplier has magnitude less than one 
up to r pa 40.3 (till 'A' in Fig. Hj), hence the limit cycle (periodic orbit) remains stable till 
r pa 40.3. The system undergoes a Niemark-S acker (NS) near r pa 40.3 as illustrated in 
Fig. [5(a). As a result of NS bifurcation, a second frequency (incommensurate with the first 
one) is generated, and the phase space trajectories show a transition from periodic orbits to 
quasiperiodic orbits. The quasiperiodic state is shown as a green patch labeled 'QP' in the 
bifurcation diagram (Fig. [3]). The phase space trajectories in this regime lie on a torus as 
illustrated in Fig. E](a) for r = 42. The power spectral density of the mode |Wioi| shown in 
Fig. [3(b) have two leading frequencies whose frequency ratio is is approximately 3.099. To 
identify the leading frequency, the power spectral density for the periodic state is shown in 
Fig-EKa). 

As the bifurcation parameter r is increased further, the two frequencies get phase locked 
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FIG. 4: The largest Floquet multiplier of the low-dimensional system as a function of r. The 
largest Floquet multiplier is greater than 1 beyond r = 40.3 (point 'A') except in the 'BC window 
in which we observe periodic states. 

to yield a periodic solution. This occurs in a narrow window in Fig. [3] and its lower inset. 
The phase space trajectories still lie on a torus but they form a periodic orbit as shown 
in Fig. MJo) for r = 44. The time-period of this periodic orbit is approximately 150 non- 
dimensional time units, and the ratio of the two leading frequencies /1//2 is approximately 
3.152 (~ 167/53). In Fig.[5]^b) we show the movement of the Floquet multipliers in the phase 
locked region. During this movement of the largest Floquet multipliers the two frequencies 
get phase locked. The power spectrum of the mode W101 for the phase-locked state is 



10] and Curry et. al. [2oJ| report /1//2 to be 



illustrated in Fig. [3(c). Gollub and Benson 
approximately 7/3 and 10/3 respectively for their phase-locked regime. Our /1//2 ~ 167/53 
is in general agreement with these earlier results. 

With a further increase in the bifurcation parameter, the system becomes chaotic at 
around r = 44.2. The route to chaos is similar to that in Curry- Yorke model (§VIII.3 



3l| ) where chaos appears after quasiperiodicity in T (2-Torus) and phase locking. In the 



bifurcation diagram (Fig. [3]), the chaotic region is shown by colored pink patch labeled 
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FIG. 5: Representation of the Floquet multiplier (FM) for various scenarios. Blue points indicate 
initial stage, red points indicate intermediate stage, and grey points indicate the final stage of 
the movement of the FM. (a) The FM crosses the unit circle through NS bifurcation creating a 
quasiperiodic state, (b) The motion of the FM during the phase-locked regime. Note that the 
largest FM remain outside the unit circle in this regime, (c) The largest FM moves into the unit 
circle resulting in a periodic solution, (d) The largest FM crosses the unit circle again creating a 
quasiperiodic solution. 

as 'CH'. As shown in Fig. [7(d) the power spectrum of the mode W W i is broad indicating 
chaotic nature of the attractor. At around r = 45.9, the size of the chaotic attractor suddenly 
increases as a result of an 'interior crisis'. This feature is illustrated in Fig. [8] where we plot 
the phase space projection on the |Wioi| — l^icd plane at r = 45.8 and r = 45.9. 

The chaotic state described above exists till r < 46.2 after which we observe periodic 
solutions (the red curves in Fig. [3]) that emerge from an inverse NS bifurcation. This inverse 
NS bifurcation is illustrated in Fig. 0(c) wherein the largest Floquet multipliers enters into 
the unit circle making the limit cycles stable. This stable time-periodic orbit continues till 
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(a) 
r = 



42 



(b) 

r = 44 




FIG. 6: A three-dimensional phase space projection of the phase space trajectories, (a) At r = 42 
the phase space trajectories fill the torus (quasiperiodic). (b) At r = 44 the system is in a phase- 
locked state, and the phase space trajectory does not fill the torus (limit cycle). 

r = 47.4 at which point another pair of Floquet multiplier again crosses the unit circle in a 
forward NS bifurcation (see Fig. 0(d)) giving rise to a quasi-periodic state. In the periodic 
window the largest Floquet multiplier is 1 as illustrated by the 'BC window in Fig. [H This 
quasi-periodic state subsequently becomes chaotic at r = 48.4 that continues for higher 
values of r. A zoomed portion of this regime of r is shown in the upper inset of Fig. [3j We 
also note that the size of the chaotic attractor is much larger than the QP attractor. This 
feature can be explained using 'attractor-merging crisis' to be discussed below. 

We discussed earlier that in the band r = 46.2 — 48.4 the low-dimensional model has a 
periodic and a quasiperiodic attractor. However, in the same band of r, a different set of 
initial conditions yield another attractor which is chaotic (see Fig. [TO]) . These two attractors 
have been shown in Fig. E^a), with the green region as the QP attractor and the pink region 
as the chaotic attractor. At r = 48.4 these two attractors merge through 'attractor-merging 
crisis' and form a single large attractor shown in Fig.[9^b). The size of the resulting attractor 
is much larger that the original QP attractor but similar to that of the chaotic attractor of 
Fig. [TOJ 

To ascertain the chaotic nature of the solutions obtained in our low-dimensional system, 
we compute the Lyapunov exponents associated with the various solutions presented in 



14 



, x 10" 



(a) r = 35 



0.0001 
10 



q^Q-eO.OS 0.1 0.15 



0.2 



0.25 



0.3 



0.35 



m 
c 

0> 
Q 

£ 0.0001 1 

o 

CD 
Q. 

w 

<D 
g 

o 

Q- 

0.0001 
10 



f2 



0.4 0.45 
(b) r = 42 



10 



0.01 





Jill 



(d) r = 45 



0.5 



_ 6 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
x10 



I 


I I I 
f2~ fix 53/167 


I 

f1 


I I 


1 1 

(c) r = 44 


I.J 


V - I I A I 


I 


1 





x10 -s0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

Frequency 



FIG. 7: The power spectral density of the mode W%oi for various dynamical states: (a) periodic, 
(b) quasiperiodic (/1//2 ~ 3.099), (c) phase locked (/1//2 ~ 3.152 ~ 167/53), (d) chaos. 

Figs. [3] and [10J The three largest Lyapunov exponents of our system corresponding to the 
attractors in Figs. [3] and [10] are shown in Figs. (TT] and [12] respectively. There is at least 
one zero Lyapunov exponent throughout the range consistent with the fact that our system 
is autonomous. In the chaotic regions described earlier there are two positive Lyapunov 
exponents clearly distinct from the zero exponent ascertaining the chaotic nature of the 
solutions. The largest Lyapunov exponents in Fig. [11] is zero for r = 46.2 — 48.4 that 
corresponds to the periodic and quasiperiodic window shown in Fig. [3] 

IV. DISCUSSIONS AND CONCLUSIONS 



In this paper we present a bifurcation analysis of a 30-mode model for the Prandtl number 
P = 6.8 (a typical large-Prandtl number fluid) and aspect ration 2V2. In our bifurcation 
analysis we observe various patterns: steady rolls (r = 1 — 27.6), time-periodic rolls (r = 
27.6 — 40.3), quasiperiodicity (r = 40.3 — 43.8), phase locking (r = 43.8 — 44.2), and chaos 
(r > 44.2). The route to chaos is similar to that of Curry- Yorke model 3lJ where chaos 
occurs after quasiperiodicity and phase locking. Periodic and quasiperiodic rolls reappear 
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FIG. 8: The phase space projection on the | Wioi | — 1^1031 plane of the chaotic attractor at (a) 
r = 45.8 and (b) 45.9. The chaotic attractor shows a sudden increase in size at r = 45.9. 

after the chaotic state in the range of r = 46.2 — 47.4 and r = 47.4 — 48.4 respectively. 
After the second quasiperiodic window the system becomes chaotic again through 'crisis'. A 
distinct feature of our low- dimensional model is that we track the fixed points, limit cycles, 
and chaotic attractors, thus getting a detailed bifurcation picture for the range of r under 
investigation. 

The above features of our 30-mode model closely resemble some of the past experimental 
results on large-Prandtl number convection namely that of Gollub and Benson lOj who 
observed chaos in water for various Prandtl numbers and aspect ratios. The route to chaos 
in Gollub and Benson's experiment for P = 5 and aspect ratio of 3.5 is through quasiperi- 
odicity and pha se locking. In direct numerical simulation of 3D RBC, Curry et al. {^(J 
and Yahata [2lJ observed similar transition to chaos for P = 10 and P = 5 respectively. 
Our low- dimensional model follow the same route to chaos. The range of Rayleigh numbers 
for our low- dimensional model is quite close to the Gollub and Benson's experiments and 
Curry et a/.'s and Yahata's DNS. Thus our low- dimensional model appears to capture the 
dynamics of 3D RBC responsible for transition to chaos through quasiperiodicity and phase 
locking. 
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FIG. 9: The phase space projection on the |Wioi| — | W103 1 plane at (a) r = 48.3 and (b) 48.5. (a) 
At r = 48.3, two attractors, a chaotic attractor (pink patch) and a quasi-periodic attractor (green 
patch) coexist. These two attractors are accessible through two different initial conditions, (b) At 
r = 48.5 the two attractor merges to generate a larger chaotic attractor. 

Our low-dimensional model also exhibits coexistence of several attractors. In the window 
of r = 46.2 — 48.4, the system has periodic and quasiperiodic attractor along with a chaotic 



attractor. Coexistence of patterns and different attractors have been observed earlier 



32] 



Another novel feature of our low- dimensional model is that it reproduces reap pearance of 



20| and Paul et 



periodic rolls after chaos, a feature observed in the 2D DNS of Curry et al. 
al. |26[], albeit at a much different r value. This feature appears to be due to the delay of 



secondary instabilities in 2D DNS compared to 3D DNS. 



Gollub and Benson 



101 ] also reported chaos in their large-Prandl number RBC through 
period-doubling, generation of three frequency (quasiperiodicity), and intermittency for dif- 
ferent sets of Prandtl numbers and aspect ratios. Our preliminary investigation for P = 10 
and aspect ratio of 2y/2 appears to indicate intermittency, however, we need to study this 
phenomena more carefully. Some of the features reported by Gollub and Benson could pos- 
sibly be captured by our model by varying the aspect ratio, a topic to be investigated in 
future. Further work is required in construction and analysis of more refined models for 
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FIG. 10: Bifurcation diagram for the low-dimensional model with initial condition different from the 
one used for generating Fig. [3l 'FP' stands for fixed point state (blue), 'OS' stands for oscillatory 
state (red), 'QP' stands for quasi-periodic state (green) and 'CH' stands for chaotic state (pink). 
'NS' stands for the Neimark-S acker bifurcation point. The present attractor and that of Fig. [3] 
differ for r = 46.2 - 48.4. 

RBC. 
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FIG. 11: Three largest Lyapunov exponents of the low-dimensional model corresponding to the 
attractors in Fig. [3J 
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FIG. 12: Three largest Lyapunov exponents of the low-dimensional model corresponding to the 
attractors in Fig. [10j 
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APPENDIX A: NUMERICAL SEARCH FOR PERIODIC SOLUTIONS OF UN- 
KNOWN PERIODS THROUGH FIXED POINTS OF A MAP 



Consider a system of first order autonomous ODEs given by 

£+ f{a,x) = 0, (Af) 

where, x G M n , a is a parameter and f(a,x) is a known vector valued function and is such 
that Eq. (IA1I) has a periodic solution. For each value of the parameter a, we seek the initial 
conditions A corresponding to a periodic solution and the time period T of the periodic 
solution. Hence, if we numerically integrate the above equation, i.e., Eq. IA1|) with x(0) = A 
till time T, we should have 

x(T) -1=0. (A2) 



Equation (1A2|) gives us n algebraic equations. However, we have n + 1 unknowns, viz. the n 
components of A and the time-period of the periodic solution T, and hence, we require one 
more equation. This equation is obtained by putting a constraint that the initial condition 
A correspond to the extremum of one of the components, e.g., x\. Accordingly, we will 
require 

Xl (T) = h(a,x(T)) = 0. (A3) 

Denoting 

y = 

we can write Eqs. ( IA2I) and (1A3I) as 

9(y) = 0. (A4) 

Equation ( 1A4I) can now be solved numerically using the Newton-Raphson method. Note 
that evaluation of g in the solution procedure involves background numerical solution of the 
ODEs (Eq. (1A1I) ). This numerical integration can be avoided if an analytical solution were 
available for the ODEs such that the map g would be known analytically. However, the 
basic principle of constructing the map y = g{y) whose fixed points gives us the appropriate 
initial conditions and the time-period of the unknown periodic solution of the ODEs remains 
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the same. 
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